{
  "metadata": {
    "kernelspec": {
      "name": "xpython",
      "display_name": "Python 3.13 (XPython)",
      "language": "python"
    },
    "language_info": {
      "name": ""
    }
  },
  "nbformat_minor": 5,
  "nbformat": 4,
  "cells": [
    {
      "id": "02898a7f-46a3-4742-85b5-efd3f368dd60",
      "cell_type": "code",
      "source": "import numpy as np\nfrom qutip import Qobj, sigmax, sigmay, sigmaz, mesolve, Options\n\n# Define physical constants and conversion factors\nh = 6.626e-34  # Planck's constant (J·s)\nc = 3.00e8     # Speed of light (m/s)\nhbar = h / (2 * np.pi)\n\n# Parameter selection based on LLLT wavelengths\nwavelengths = np.array([630e-9, 660e-9, 810e-9, 980e-9])  # Common LLLT wavelengths (m)\nomega_0 = 2 * np.pi * c / wavelengths  # Transition frequency (Hz)\n\n# Scale Rabi frequency (Ω) based on fluence/irradiance\nfluence = 10.0  # Typical fluence ~10 J/cm²\nirradiance = 0.1  # Irradiance ~100 mW/cm²\narea = 1e-4      # Beam area ~1 cm²\npower = irradiance * area\nenergy_per_photon = h * c / wavelengths\nphoton_flux = power / energy_per_photon\nOmega = 2 * np.pi * photon_flux / 1e6  # Rabi frequency (MHz), scaled to photon flux\n\n# Biological decoherence rate (γ) based on physiological conditions\ngamma = np.array([0.05e6, 0.1e6, 0.2e6])  # Decoherence rates (MHz), reflecting pH, temperature effects\n\n# Define TLS Hamiltonian\nsigma_x, sigma_y, sigma_z = sigmax(), sigmay(), sigmaz()\nH0 = 0.5 * omega_0[0] * sigma_z  # Base Hamiltonian with first wavelength\nH1 = 0.5 * Omega[0] * sigma_x    # Interaction Hamiltonian\n\n# Time parameters\ntlist = np.linspace(0, 50e-6, 500)  # Time array (µs)\n\n# Lindblad master equation setup\ndef lindblad_deriv(rho, t, args):\n    H = H0 + args['Omega'] * H1\n    return -1j * (H * rho - rho * H.dag()) + sum(g * (s * rho * s.dag() - 0.5 * (s.dag() * s * rho + rho * s.dag() * s)) for g, s in zip(args['gamma'], [sigma_x, sigma_y, sigma_z]))\n\n# Initial state (ground state)\npsi0 = Qobj([[1], [0]])\n\n# Solve for dynamics\noptions = Options(store_states=True)\nresult = mesolve(H0, psi0, tlist, c_ops=[np.sqrt(g) * s for g, s in zip(gamma[0], [sigma_x, sigma_y, sigma_z])], e_ops=[sigmaz()], args={'Omega': Omega[0], 'gamma': gamma[0]})\n\n# Output excited population\nexcited_pop = 0.5 * (1 - result.expect[0])",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "id": "912350cb-71c9-40c6-a940-1411f906bf80",
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    }
  ]
}